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We measure the recoil velocity as a function of spin for equal-mass, highly-spinning black-hole 
binaries, with spins in the orbital plane, equal in magnitude and opposite in direction. We confirm 
that the leading-order effect is linear in the spin and the cosine of angle between the spin direction 
and the infall direction at merger. We find higher-order corrections that are proportional to the odd 
powers in both the spin and cosine of this angle. Taking these corrections into account, we predict 
that the maximum recoil will be 3680 ± 130km s~ . 
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I. INTRODUCTION 

The field of Numerical Relativity (NR) has progressed 
at a remarkable pace since the breakthroughs of 2005 [T]- 
|3] with the first successful fully non-linear dynamical nu- 
merical simulation of the inspiral, merger, and ringdown 
of an orbiting black-hole binary (BHB) system. BHB 
physics has rapidly matured into a critical tool for gravi- 
tational wave (GW) data analysis and astrophysics. Re- 
cent developments include: studies of the orbital dynam- 
ics of spinning BHBs [4TfTT] , calculations of recoil veloci- 
ties from the merger of unequal mass BHBs [T2TH4] , and 
very large recoils acquired by the remnant of the merger 
of two spinning BHs 8, 15 -29 , empirical models re- 
lating the final mass and spin of the remnant with the 
spins of the individual BHs [30TI37] , and comparisons of 
waveforms and orbital dynamics of BHB inspirals with 
post-Newtonian (PN) predictions [3"8144"5] . 

The surprising discovery jTS], [TB] that the merger of 
binary black holes can produce recoil velocities up to 
4000 kms -1 , and hence allow the remnant to escape from 
major galaxies, led to numerous theoretical and observa- 
tional efforts to find traces of this phenomenon. Several 
studies made predictions of specific observational features 
of recoiling supermassive black holes in the cores of galax- 
ies in the electromagnetic spectrum [4"6"H5"2"] from infrared 
[53] to X-rays [SIHSfj] and morphological aspects of the 
galaxy cores [57H59] . Notably, there began to appear ob- 
servations indicating the possibility of detection of such 
effects [50Tfrj2"] , and although alternative explanations are 
possible 63 66J , there is still the exciting possibility that 
these observations can lead to the first confirmation of a 
prediction of General Relativity in the highly-dynamical, 
strong-field regime. 

Numerical simulations of the BHB problem have sam- 
pled the parameter space of the binary for different val- 
ues of the binary's mass ratio q and arbitrary orienta- 
tions of the individual spins of the holes. Two astrophys- 
ically important regions of this parameter space remain 
challenging to describe accurately by numerical simula- 
tions: the small q limit, although recent development 



of the numerical techniques have produced a successful 
simulation of the last few orbits before the merger of 
a q = 1/100 binary [67], and the near maximal spin 
limit. The most recent simulations of highly-spinning 
BHBs was of non-precessing binaries with intrinsic spins 
a = 0.95 [55]. Since BHBs with a = 1 are still elusive 
to full numerical simulations, and the configuration that 
maximizes the gravitational recoil is one that starts with 
maximally spinning BHs, with opposite spins lying on the 
orbital plane [T5] [H] , we will model these configurations 
for different values of the intrinsic spin parameter up to 
a = 0.92 (which is achievable with current techniques 
to solve initial "puncture" data) and then extrapolate to 
a = 1 using an improves version of our original empirical 

formula nana [HO]. 

In Ref. |30) we extended our original empirical formula 
for the recoil velocity imparted to the remnant of a BHB 
merger [151 116) to include next-to-leading-order correc- 
tions, still linear in the spins. The extended formula has 
the form: 
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where 77 = q/(l + q) 2 , with q — mx/m2 the mass ra- 
tio of the smaller to larger mass hole, di = Si/m 2 , the 
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index _L and || refer to perpendicular and parallel to the 
orbital angular momentum respectively, e±, £2 are orthog- 
onal unit vectors in the orbital plane, and £ measures the 
angle between the unequal mass and spin contribution to 
the recoil velocity in the orbital plane, from newly avail- 
able runs. The angle 9 is defined as the angle between 
the in-plane component of A = M(Sijm2 — Si/m\) or 
S = S\ + S% and a fiducial direction at merger (see 
Ref. [TH] technique). Phases O and ©i depend on the 
initial separation of the holes for quasicircular orbits (as- 
trophysically realistic evolutions of comparable masses 
black holes lead to nearly zero eccentricity mergers). 

The empirical formula above was obtained by as- 
suming the post-Newtonian dependence on the spin and 
mass ratio of instantaneous radiated linear momenta [69] 
where the coefficients are to be fitted by full numerical 
simulations. Second order corrections in the spin have 
been obtained recently [70] and could be added to the 
empirical formula. Here, in this paper, we will consider 
instead the particular family of configurations that lead 
to the maximum recoil [15H17L 171) . where q = 1 and 
the two spins are in the orbital plane, equal in magni- 
tude, and opposite in direction. These configurations are 
7T— symmetric, i.e. rotating the system by 180 degrees 
around the symmetry axis lead to the same configura- 
tion. This implies in particular, that only odd powers of 
the spin and the cos are involved. We will then perform 
a series of simulations that vary both the magnitude of 
the (intrinsic) spin in the range a = 0.2 — 0.92 and the 
initial angle of the individual black-hole spin and orbital 
linear momentum. 

For a first exploration of the extended spin dependence 
we consider cubic and possible fifth-order corrections [32] 
to the empirical formula (JlJ of the form 

V|| = (V hl a + V 1:3 a 3 ) cos(0 A - ©o) 

+ (V 3 ,ia + V 3 W + V^a 4 ) cos(30 A - 30 3 ), (2) 
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where Vx t i = 2K(1+t]Bk) ^ +q ^ , and the remaining terms 
are higher-order correction to Eq. 0. 

II. TECHNIQUES 

To compute the numerical initial data, we use the 
puncture approach [72] along with the TwoPUNC- 
tures [73] thorn. In this approach the 3-metric on the 
initial slice has the form j a f, = (ipBL + u) 4 5 a i,, where 
ipBL is the Brill-Lindquist conformal factor, S a b is the 
Euclidean metric, and u is (at least) C 2 on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
tpBL = 1+X)j=i m i/(2\r — 'Til), where n is the total num- 
ber of 'punctures', is the mass parameter of puncture 
i (mf is not the horizon mass associated with puncture 
i), and rj is the coordinate location of puncture i. For the 
initial (conformal) extrinsic curvature we take the ana- 
lytic form Kfj Y given by Bowen and York[?l]. We evolve 



these black-hole-binary data-sets using the LazEv [73] 
implementation of the moving puncture formalism [U [2] 
with the conformal factor W = Jx = exp(— 2(f)) sug- 
gested by [1] as a dynamical variable. For the runs pre- 
sented here we use centered, eighth-order finite differenc- 
ing in space [76] and an RK4 time integrator (note that 
we do not upwind the advection terms). 

We use the Carpet [77j mesh refinement driver to pro- 
vide a 'moving boxes' style mesh refinement. In this ap- 
proach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then 
moves these fine grids about the computational domain 
by following the trajectories of the two black holes. 

We use AHFinderDirect [78] to locate apparent 
horizons. We measure the magnitude of the horizon 
spin using the Isolated Horizon algorithm detailed in 79 . 
This algorithm is based on finding an approximate rota- 
tional Killing vector (i.e. an approximate rotational sym- 
metry) on the horizon tp a . Given this approximate Killing 
vector ip a , the spin magnitude is 

S [v] = f & a R b K ab )d 2 V, (3) 

m J AH 

where K ab is the extrinsic curvature of the 3D-slice, d 2 V 
is the natural volume element intrinsic to the horizon, 
and R a is the outward pointing unit vector normal to 
the horizon on the 3D-slice. We measure the direction of 
the spin by finding the coordinate line joining the poles 
of this Killing vector field using the technique introduced 
in [3J. Our algorithm for finding the poles of the Killing 
vector field has an accuracy of ~ 2° (see [3J for details) . 
Note that once we have the horizon spin, we can calculate 
the horizon mass via the Christodoulou formula 

m H = ^/m? r + S 2 /(4m? r ), (4) 

where m- m — \J Aj (167r) and A is the surface area of the 
horizon. We measure radiated energy, linear momentum, 
and angular momentum, in terms of tp^, using the for- 
mulae provided in Rcfs. [80] [81]. However, rather than 
using the full ip4, we decompose it into £ and m modes 
and solve for the radiated linear momentum, dropping 
terms with I > 5. The formulae in Refs. 80, 8T] are valid 
at r = 00. We obtain highly accurate values for these 
quantities by solving for them on spheres of finite radius 
(typically r/M = 50,60, • • • , 100), fitting the results to 
a polynomial dependence in / = 1/r, and extrapolating 
to I = EH [H2 [83]. Each quantity Q has the ra- 
dial dependence Q = Qo + IQi + 0(l 2 ), where Qo is the 
asymptotic value (the 0(1) error arises from the Oil) er- 
ror in rtpi). We perform both linear and quadratic fits 
of Q versus /, and take Qo from the quadratic fit as the 
final value with the differences between the linear and 
extrapolated Qo as a measure of the error in the extrap- 
olations. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1+log lapse and a modified Gamma-driver 
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shift condition [HUH]; and an initial lapse a(t = 0) = 
2/(1 + The lapse and shift are evolved with 

{d t -P l di)a = -2aK, (5a) 
d t p a = (3/4)f a - v (x a ,t)(3 a , (5b) 

where different functional dependences for rj(x a ,t) have 
been proposed in [751 I85H89"] . For the low-spin simula- 
tions we used a constant r\ = 2, while for the a = 0.92 
simulation we used a modification of the form proposed 
in |M 



r)(x a ,t) = R Q 



y/fidiWdjW 
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where we chose Ro = 1.31 [SB]. In practice we used 
a = 2 and 6 = 2, which reduces rj by a factor of 4 at 
infinity when compared to the gauge proposed by [55] , 
improving its stability at larger radii. Other values of 
(a, b) lead to an increase of the numerical noise. Note 
that this gauge was originally proposed and used for the 
non-spinning, intermediate-mass-ratio binaries. Here we 
find that the gauge is well adapted for the highly-spinning 
equal mass case, where, after the initial burst of radiation 
passes, the measured spin is found to never drop below 
a = 0.905. Due to the differences in the spurious initial 
radiation content, as well as spin-orbit effects on the total 
mass, a near merger varied from between 0.90 to 0.93 for 
the different A09Tyyy configurations (See tables [T| and [ll| 
below). 



A. Initial Data 

We used 3PN parameters for quasicircular orbits with 
BH spins (equal in magnitude and opposite in direction) 
aligned with the linear momentum of each BH (i.e. in- 
plane spins) to obtain the momenta and spin parameters 
for the Bowen-York extrinsic curvature. We then chose 
puncture mass parameters such that the total ADM mass 
was 1M. We then rotated the spins by 30°, 90°, 130°, 210° 



and 315°, to obtain a total of 6 configurations for each 
value of the intrinsic spin a. We label the configuration 
AxxTyyy where xx corresponds to the spin of each BH 
and yyy is the initial rotation of the spin directions. We 
summarize the initial data in Table [T| 

This family of configurations has larger initial separa- 
tions than the configurations in our original studies in 
|16| . In these configurations, the BHs orbit ~ 3.5 times 
prior to merger, which allows for most of the eccentric- 
ity to be radiated away before the plunge phase (where 
most of the recoil velocity is generated). This provides 
for an accurate description of the plausible astrophysical 
maximal recoil scenario. 

TABLE I: Initial data parameters for the non-rotated config- 
urations. The initial puncture positions are ±(i, 0, 0), mo- 
menta are ±(0,p, 0), and spin ±(0, S, 0). The remaining con- 
figurations are obtained by rotating the spins, keeping all 
other parameters the same. 



Config 


X 


P 


s 


m p 


A02T000 


3.878113 


0.117404 


0.051314 


0.479782 


A04T000 


3.879566 


0.117405 


0.102627 


0.454076 


A06T000 


3.881979 


0.117407 


0.153936 


0.403550 


A08T000 


3.885342 


0.117409 


0.205241 


0.301026 


A09T000 


3.887375 


0.117411 


0.230891 


0.172120 



III. RESULTS AND ANALYSIS 

In order to analyze our results for different initial ori- 
entations of the spin that span the 0— dependence, we 
use the techniques detailed in [IS]. For each a we fit 
the results of the recoil as a function of angle to form 
Kocoii = V\ cos(d - 0i) + V 3 cos[3(0 - 3 )], where is de- 
fined to be the angle of the spin direction (of the first BH) 
near merger (at a fiducial radial separation of r = 1.2) 
and the spin direction of the corresponding AxxTOOO con- 
figuration (we cannot simply use the initial spin direction 
differences because spin-orbit effects for larger spins make 
this approximation inaccurate). The radiated energy and 
recoil from each simulation is given in Table |Hj 



We then fit V\ and V3 to the functional forms V\ = 
Vi^a+Vi^a 3 and V3 = V^ia+V^aO! 3 . A summary of the 
fits is given in Table III Note that Vi 1 is related to the 
parameter K in our empirical formula ^ by K — 16V\^. 
Here we find K = 58912 ±43, where the error is obtained 
from the fit and likely underestimated the true error in 
this quantity. Previously we found K — (6.0 ±0.1) x 10 4 , 
which agrees reasonably well with the new value [IB] US] • 
We also include fits where the linear term in V3 and the 
cubic term in V\ are set to zero, as well as a fit of V3 
to V 3 . 3 a 3 + V 3}5 a 5 . We note that a cubic term in V\ is 
expected since cos 3 = 3/4 cos + 1/3 cos 30, and hence 
cubic corrections of the for a 3 cos 3 will contribute to 



the cos0 dependence. On the other hand, a linear de- 
pendence in cos 30 is not expected. 

The form of the fitting above was first proposed in 
[32] as a generic expansion, where it was applied to data 
sets with constant a. Here we compare results from five 
different values of the intrinsic spin in the range a = 0.2— 
0.92, to obtain an accurate model of the a dependence. 



In Figs. [T][5] we show the angular fits for each set of 
configurations. Note that the spin-orbit coupling effects 
are strongest for the A09Tyyy configurations, as is ap- 
parent by the relative translation of two configurations 
towards the same final angle. 
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TABLE II: The radiated energy, recoil velocity, and angle between the spins for the AxxTyyy configuration at merger and the 
corresponding AxxTOOO configuration, AO. Note the substantial rotations apparent in the A09Tyyy configurations due to spin 
orbit interactions. 



Config 


SE 


Vrecoil 


A6 


SJ Z 


A02T000 


0.03583 ± 0.00020 


551.95 ±0.83 





0.2727 ±0.0032 


A02T030 


0.03575 ± 0.00020 


225.49 ± 0.86 


30.7 


0.2724 ±0.0031 


A02T090 


0.03562 ±0.00019 


-482.10 ±0.13 


88.6 


0.2721 ± 0.0031 


A02Tf30 


0.03574 ± 0.00019 


I AO 1 A O A 

—721.03 ± 0.34 


127.3 


A A^TA^T 1 A A A O ~1 

0.2727 ± 0.0031 


A02T210 


0.03573 ± 0.00020 


-234.14 ±0.76 


210.0 


0.2723 ± 0.0032 


A02T315 


0.03577 ±0.00019 


730.10 ±0.32 


312.4 


0.2729 ± 0.0030 


A04T000 


0.03665 ± 0.00022 


1200.79 ±2.21 





0.2768 ± 0.0031 


A04T030 


0.03625 ± 0.00022 


529.09 ±2.12 


33.7 


0.2747 ±0.0031 


A04T090 


0.03574 ±0.00020 


-764.08 ± 0.39 


85.0 


0.2740 ± 0.0028 


A a A rri i oa 

A04Tf30 


0.03620 ± 0.00021 


■I OAA AA 1 A AA 

— 1390.22 ± 0.89 


126.0 


A ATTA 1 A AAAA 

0.2759 ± 0.0029 


A04T210 


0.03633 ±0.00021 


-637.049 ± 1.97 


209.2 


0.2755 ± 0.0030 


A04T315 


0.03628 ± 0.00021 


1424.17 ± 1.02 


311.2 


0.2762 ± 0.0029 


A06T000 


0.03803 ± 0.00027 


2001.06 ±3.5 





0.2834 ±0.0032 


A06T030 


0.03740 ± 0.00026 


1087.24 ±4.8 


36.2 


0.2798 ± 0.0032 


A06T090 


0.03595 ± 0.00025 


-870.93 ± 2.4 


87.3 


0.2758 ± 0.0026 


A06Tf30 


0.03669 ± 0.00026 


— 1944.04 ± 1.5 


127.7 


0.2791 ± 0.0029 


A06T2fO 


0.03751 ± 0.00026 


-1212.62 ± 4.4 


212.4 


0.2807 ±0.0031 


A06T315 


0.03679 ± 0.00026 


1984.74 ± 1.3 


310.3 


0.2797 ±0.0029 


A08T000 


0.03996 ± 0.00039 


2651.75 ±7.54 





0.2912 ±0.0027 


A08T030 


0.03941 ± 0.00037 


1917.03 ±8.29 


24.3 


0.2888 ± 0.0027 


A08T090 


0.03677 ±0.00034 


-445.83 ± 1.10 


70.9 


0.2791 ± 0.0028 


A08T130 


0.03733 ± 0.00038 


-2412.2 ±4.03 


119.5 


0.2823 ± 0.0027 


A08T210 


0.03941 ± 0.00036 


-1919.86 ±8.26 


204.3 


0.2887 ±0.0027 


A08T315 


0.03771 ± 0.00038 


2568.87 ±4.49 


306.3 


0.2838 ± 0.0027 


A09T000 


0.04026 ± 0.00057 


89.74 ± 1.00 





0.3063 ± 0.0042 


A09T030 


0.04143 ± 0.00055 


3240.96 ± 17.34 


101.5 


0.3011 ±0.0051 


A09T090 


0.04062 ± 0.00057 


1859.42 ± 15.47 


147.4 


0.2951 ± 0.0070 


A09T130 


0.03784 ± 0.00054 


-759.16 ±0.85 


190.6 


0.2846 ± 0.0058 


A09T210 


0.04144 ±0.00055 


-3239.25 ± 17.24 


281.7 


0.3012 ± 0.0050 


A09T315 


0.03917 ±0.00056 


-205.59 ±2.92 


355.3 


0.2919 ± 0.0067 



In Table III we provide the fitting constants Vij as- 
suming the spin of the A09Tyyy was 0.92. in actuality, 
the spin varied between configurations. In Table |IV| we 
provide fitting parameters for Vij if we take the value of 
a for these configurations to be a = 0.9 (the expected 
value when neglecting effects due to the initial radiation 
content), a = 0.91, and a = 0.92 (which approximates 
the average value of a over all configurations). We find 
that setting a — 0.92 gives the best fit for the dominant 
Vj.,1 term. However, we note that these fits do indicate 
that the nonleading term and V^i term may be zero. 
We therefore also provide fits assuming these two terms 
vanish. Fits to V\ strongly prefer a = 0.92 over the 
smaller values. We note that the sign of Vi.,3 changes if 
we assume smaller values of a for the A09Tyyy configu- 
rations. 



FIG. 1: Fit of the recoil versus angle for the a — 0.2 configu- 
rations. 
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While arguments based on post-Newtonian scaling do 
not seem to indicate the presence of an a cos 30 term, 
our results indicate that this term is present. This may 
indicate an error in V3 for a = 0.2. If we exclude this 
data point, then we can fit to reasonably well to either 
V^ia ± V3 j3 a 3 or V 3 _ 3 a 3 ± V^a 5 . Further exploration in 
the small a regime is required. In Figs [6] we compare fits 



of V3 and find that the best fit is to V3 = aV^i ± a 3 V3,3- 
On the other hand, as seen in Fig.[7J there is no significant 
difference apparent in the fits of V\ to Vi^a + Vi^a 3 and 
Vi = Vi ia. 
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TABLE III: Fits of the recoil to the functional form Wkick = 
Vi cos(# — 9i) ± V3 cos[3(# — S3)]. Note that angles are mea- 
sured in degrees. The reported errors come from the nonlinear 
least-squares fit of the data and are underestimates of the ac- 
tual errors. Note the average value of a for the A09Tyyy 
configuration near merger was a ~ 0.92, which was the value 
used to obtain the fits to Vij. 



a 


Vi 


9i 








>fc zfc ^ 


0.2 


737.70 ± 0.12 


221.8002 ± 0.0010 


0.4 


1472.59 ± 0.06 


215.6909 ± 0.0011 


0.6 


2204.98 ± 0.56 


205.117 ± 0.015 


0.8 


2935.93 ±0.65 


206.658 ±0.013 


0.9 


3376.3 ± 7.5 


91.02 ±0.11 


a 


v 3 


03 





**** 




0.2 


4.23 ±0.12 


279.62 ±0.65 


0.4 


12.0838 ± 0.024 


37.790 ± 0.049 


0.6 


31.63 ±0.55 


152.72 ±0.38 


0.8 


69.21 ±0.74 


38.01 ±0.22 


0.9 


95.5 ±2.4 


36.7 ± 1.5 


Vl,i 


3681.77 ±2.66 




Vi, 3 


-15.46 ±3.97 




V 3 ,i 


15.65 ±3.01 




v 3 , 3 


105.90 ±4.50 





TABLE IV: Fits Vi and V 3 to the form Vx = Vi,ia + Vi, 3 a 3 
and V 3 = V 3 ,ia + V3,3a 3 , as well as V 3 = V 3l3 a 3 + V 3 ^a 5 . 
For the A09Tyyy configurations we take a = 0.9, 0.91, and 
0.92, which accounts for the expected value of a for these 
configuration, the actual average value observed, and a spin 
between these two values, as explained in the text. 5 2 is the 
average of the square of the error in the fit. Note that fits to 
the dominant Vi term strongly prefer a — 0.92 over smaller 
values, while fits to the subleading V 3 prefer a = 0.9. 



a (A09Tyyy) 


Vi,i 


Vl,3 


8 2 


0.92 


3681.77 ±2.66 


-15.46 ±3.966 


1.21 


0.91 


3658.21 ± 20.74 


49.16,31.47 


71.13 


0.90 


3634.85 ± 41.09 


115.31,63.40 


270.25 


a (A09Tyyy) 


Vz,i 


v 3 , 3 


S 2 


0.92 


15.65 ±3.01 


105.90 ±4.50 


1.55 


0.91 


13.68 ± 1.82 


111.15 ±2.77 


0.55 


0.90 


11.75 ± 1.14 


116.45 ± 1.77 


0.21 


a (A09Tyyy) 


Vi,i 




5 2 


0.92 


3672.08 ± 1.84 





5.80 


0.91 


3688.56 ±8.23 





114.53 


0.90 


3704.98 ± 17.17 





493.78 


a (A09Tyyy) 




v 3 , 3 


S 2 


0.92 





127.74 ± 3.96 


12.02 


0.91 





130.60 ±3.36 


8.29 


0.90 





133.45 ±2.90 


5.73 


a (A09Tyyy) 


v 3}3 


Vs.s 


S 2 


0.92 


172.55 ± 10.20 


-58.98 ± 13.21 


2.01 


0.91 


167.45 ± 10.87 


-49.54 ± 14.38 


2.09 


0.90 


161.80 ± 12.15 


-38.88 ± 16.43 


2.389 



FIG. 2: Fit of the recoil versus angle for the a — 0.4 configu- 
rations. 
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FIG. 3: Fit of the recoil versus angle for the a — 0.6 configu- 
rations. 
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IV. CONCLUSION 

Using the enhanced recoil formula for the "maximum 
kick" configurations, we predict that the maximum re- 
coil will be 3680 ± 130km s _1 , where the error in the 
prediction is due to the possibility of the higher-order 
effects producing recoils in the same direction or oppo- 
site direction of the dominant linear contribution. We 



FIG. 4: Fit of the recoil versus angle for the a — 0.8 configu- 
rations. 
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FIG. 5: Fit of the recoil versus angle for the a — 0.92 config- 
urations. 
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FIG. 6: A comparison of fits of V3 to V3 = Q.V3 1 + a 3 V3,3 
(solid), Vz = a 3 V 3 , 3 + a 5 V3,s (dotted), V 3 = a 3 V 3 , S (dot- 
dashed). The first fit is the best. In all cases the spins for the 
A09Tyyy configurations were assumed to be a = 0.92 

V3 
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also established a model for higher-order dependences 
on the spin in the recoil formula. These results are 
particularly relevant for the interpretation of observa- 
tions of emission lines in AGNs displaying displace- 
ments between narrow and wide emission lines of the 
order of thousands of kilometers per second. In partic- 
ular in Ref. [HI] a 1200 km/s offset velocity was mea- 
sured (CXOCJ100043.1+020637). A 2650 km/s recoil- 



ing supermassive black hole could explain the observa- 
tions (SDSS J092712+294344) of Ref. [60]. While in 
Ref. [BJ (SDSS J105041+345631) and in Ref. [92] (SDSS 
J153636+044127) there is speculation that 3500 km/s re- 
coiling black holes are responsible for these features in the 
spectra. While none of those cases effectively surpasses 
the maximum recoil velocity determined here, they came 
close enough for the probability of actually observing this 
event to be very low [53] thus leading to the question 
about what are the astrophysical mechanisms responsi- 
ble of generating such large differential velocities [SH [23 • 

FIG. 7: A comparison of fits of Vi to Vi = aVii + a 3 Vi 3 
(solid), V 3 = a Vi,i (dotted). There is no significant differ- 
ences between the fits. In all cases the spins for the A09Tyyy 
configurations were assumed to be a — 0.92 
vi 
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